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ABSTRACT 

We explore structure formation in the dark ages (z ~ 30 — 6) using two well-known methods for 
initializing cosmological TV-body simulations. Overall, both the Zel'dovich approximation (ZA) and 
second order Lagrangian perturbation theory (2LPT) are known to produce accurate present-day dark 
matter halo mass functions. However, since the 2LPT method drives more rapid evolution of dense 
regions, it increases the occurrence of rare massive objects - an effect that is most pronounced at high 
redshift. We find that 2LPT produces more halos that could harbor Population III stars and their 
black hole remnants, and they produce them earlier. Although the differences between the 2LPT and 
ZA mass functions are nearly erased by z = 6, this small boost to the number and mass of black 
holes more than doubles the reionized volume of the early Universe. We discuss the implications for 
reionization and massive black hole growth. 

Subject headings: cosmology: theory — galaxies: high-redshift — intergalactic medium — physical 
data and processes: black hole physics — physical data and processes: radiative 
transfer — methods: numerical 



1. INTRODUCTION 

The first generation of stars, so-called Population III 
stars (Pop III), form from pristine gas within dark matter 
halos (DMH) at very high r edshifts (jCouchman fc Reesl 
119861 : iTegmark et~atl 119971 ). Both one-zone calcula- 
tions and cosmological simulations suggest a top-heavy 
Pop HI IMF with an uncertain characteristic mass 
(e.g. Bromm et~aDl999l lOmukaill2QQ(l lAbel et allkooa 
ISmith et all l2009h . The uncertainty hinges on recent 
simulations that showed the fragmentation and multi- 
plicity may occur at p r otostellar densitie s (jTurk et al.l 
[20091 : IClark et al.l l201ll IGreif et alJl2012h . and on the 
idea that protostellar radiative feedback may limit mass 
accre tion to ~ 40 M (Hosok awa et al.ll2011l : iStacv et al.l 
120121 ). Massive Pop III stellar lifetimes are ~ 2-20 Myr, 
after which they either enrich the intergala ctic medium 
(IGM ) with metals through supernovae (|Heger et al.l 
120031 ) or collapse directly into black holes (BH) that 
could grow into the sup ermassive ones at galac- 
tic centers jMadau & Rees 2001; Volonte rTeTall 120031: 



„ Tanaka fc H aiman 2 varez et all 

20091 : Ueon et al.ll2011b iMicic et al.ll201lD . While the Pop 
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III IMF is still uncertain, it is likely that DMHs do host 
Pop III stars, and a significant fraction produce seed BHs. 

To form the first generation of stars, primordial gas 
collapses within the DMH center via two separate cool- 
ing channels - mo l ecular hydrogen and atomic hydrogen. 
iTrenti fc Stiavellil (|2009[ ) found the minimum DMH mass 
to form cold dense core through molecular hydrogen cool- 
ing: 
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6.44 x lO^oJ^ 457 
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31 
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(i) 

where J21 relates to the Lyman- Werner flux, F* = 
47rJ2ilO -21 ergs -1 cm -2 . It is this flux that is respon- 
sible for photo-dissociating molecular hydrogen. To cal- 
culate F*, we adopt th e comoving Lyman- Wern er pho- 
ton density, nLw, from ITrenti fc Stiavellil (l2009h . which 
assumes stellar photon yields from iSchaererl ([2003): 



^lw(^) = 7 x 10' 
and 

J21 = 1.6 x 10~ 65 



10 3.3-3.3(l+z)/ll Mpc - 




(2) 



(3) 



For J 2 i - 5 at z = 10, M min - 5 x 10 8 M . Q 

On the other hand, all DMHs with virial temperatures 
> 10 4 K trigger efficient atomic hydrogen cooling, in- 
dependent of the Lyman- Werner background. A virial 
temperature of 10 4 K translates into a DMH mass of: 



M vir = 7.75 x 10 b Mr7 



1 



31 



-1.5 



(4) 



Regardless of the cooling mechanism, the DMHs that 
host the first generation of stars are by far the most mas- 
sive structures at this early epoch. For example, a DMH 
of ~ 10 8 M is a ~ 5 — a peak at z = 20. These rare 
and massive peaks are the first collapsing structures in 
the Universe, and are the most likely to be plagued by 

4 Note that neutral HI in the IGM can reduce th e effective LW 
flux from equation [2] by an order of magnitude ( Haim an et al.l 
2000); this should not affect our results since we are concerned 
with the difference between two otherwise identical volumes. 
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numerical transients and initial value issues in any cos- 
mological simulation. The purpose of this paper is to 
compare the effect on the Pop III era of two different tech- 
niques for initializing a cosmological iV-body simulation: 
the Zel'dovich approximation (ZA), and second order La- 
grangian perturbation theory (2LPT) on the Pop III era. 
Below, we mer ely outline th e cosm ological initialization 
techniques; see iScoccimarrol (jl998l ) for a comprehensive 
overview. 

The particle positions in an TV-body simulation are 
generated from a density spectrum consistent with the 
Cosmic Microwave Background (CMB). Since TV-body 
simulations start at much lower redshifts (z sta rt ~ 50 — 
250), the particle evolution up to the starting redshift is 
estimated using a displacement field. In the Zel'dovich 
approximation, the displacement field assumes a linear 
evolution from recombination to z sta rt so the particle 
displacement is determined by the linear overdensity and 
linear growth factor D±. In 2LPT, the displacement field 
also includes a second order growth factor D2 ~ —ZD\jl 
and accounts for some of the non-linear density evolution. 

Since ZA is accurate only to first order, any higher 
order growing modes will be incorrect. In addition, 
in ZA the non-linear decaying modes, called transients, 
damp away with 1/a, while in 2LPT, transients damp 
much more quickly, as 1/a 2 . This implies that the col- 
lapse time of true structures in an TV-body simulation 
will be acc urate after fewer e-folding times when using 
2LPT (e.g.. iScoccimarr o 1998; Crocce et al.ll2006tlJenkinsl 
2010]). Thus, for same starting redshift, simulations ini- 
tialized with 2LPT will capture the formation of very 
high redshift halos more accurately compared to simu- 
lations initialized with ZA. In general, high-sigma peaks 
are supressed in ZA, and the effect is larger at high red- 
shift. Of course, given enough expansion factors, the 
transients eventually damp out for all but the most mas- 
sive DMHs. The standard lore is that a ZA simulation 
must only start earlier than 2LPT; indeed, by z = 1 there 
is only a ~ 10% difference in the D MH number density 
above 10 14 M ([Crocce et all 12000 ). This small differ- 
ence in the present day DMH mass function is perhaps 
why the Zel'dovich approximation has been adopted so 
widely. 

2. SIMULATIONS AND BLACK HOLE MODELING 

2.1. Simulations 

In any cosmological simulation, the goal is to displace 
the particles only slightly from their initial positions. 
The rms displacement, A rms is: 

47T f kny 

A rms = - P(k, Z start ) dk, (5) 

where kf = 27r/Lbox is the fundamental mode, Lbox * s 
the simulation box-size, k ny = N/2 x kf is the Nyquist 
frequency for an N 3 simulation, and P(k, z start) is the 
power spectrum at the starting redshift, z start- If A rms 
is larger than the mean interparticle separation A p = 
Lbox/N, then "orbit crossing" occurs and invalidates 
the accuracy of the initial conditions. Strictly speak- 
ing, imposing A rms /A p = 0.1 would be ideal in that it 
makes orbit crossing a ~ 10-sigma event. The challenge 
in small volume, high resolution simulations is that A p 



is small and the starting redshift must be very high to 
satisfy this A rms /A p constraint. For instance, to sat- 
isfy A rms /A p ~ 0.1, a 10 IrC 1 Mpc, 512 3 simulation 
would need z s tart ~ 799. However, at very high red- 
shifts the matter distribution is quite smooth, and the 
net f orce may be domi nated by numerical round-off er- 
rors (iLukic et al.ll2007[ ) - this will suppress small scale 
structure. Therefore, we adopt A rms /A p = 0.25, which 
sets Zstart = 300. We note that this criterion is rarely 
mentioned in the literature. Several small volume, high 
resolution simulations appear to have A rms /A p > 1.0, a 
clear sign that the initial conditions are not valid. Since 
canonical Pop Ill-hosting halos are expected to appear 
at z ~ 30, this z s tart allows 10 expansion factors to oc- 
cur, which reduces transients from the initial conditions 
in 2LPT by a factor of 100. 

We initialize six 512 3 , 10 h~ Y Mpc, dark matter-only 
cosmological volumes o f a ACDM Universe and evolved 
them with Gadget-2 (jSpringel et al.l l2001ri iSpringell 
[2QQ5h from 

Zstart — 300 to z — 6. The six simulations di- 
rectly contrast the initialization technique, because each 
pair samples the volume identically from the CMB trans- 
fer function, and only displaces these identical initial po- 
sitions to the same starting redshift using either 2LPT 
or ZA . Each volume is initi alized using WMAP-5 param- 
eters ([Komatsu et all 12009( 1. 

Given that we are concerned with the evolution of our 
volumes at very high redshift, care was taken to inte- 
grate the positions and velocities with high accuracy. We 
adopted many of the tree code p arameters of the Coyote 
Universe ([Heitmann et al.ll2010l ). PMGRID, which de- 
fines the Fourier grid is 1024, ASMTH, which defines the 
split between long and short-range forces, is 1.5 times 
the mesh cell size, RCUT, which controls the maximum 
radius for short-range forces, is 6.0 times the mesh cell 
size, the force accuracy is 0.002, the integration accuracy 
is 0.00125 and the softening is 1/40A P = 0.5ft- 1 kpc. 
The particle mass for all simulations is 5.3 x lO 5 ft _1 M0. 

2.2. Black hole growth and accretion luminosity 

In post-processing, we i dentified halos with a t least 20 
particles using SUBFIND (ISpringel et al.ll2001aD and con- 
structed a halo mergertree ([Sinha fc Hollev-Bockelmannl 
12012( 1 . We seeded 100 fc^M© BHs in halos with the mass 
criteria defined in the previous section^. Once a BH seed 
is sown, we do not allow another BH to form within that 
halo. The merger tree now allows us to track the growth 
and merger of BHs as well. We assumed a simple BH 
growth prescription: Eddington-limited accretion trig- 
gered by a major merger, where the BHs merge promptly 
after D MHs merge, and the re is no gravitational wave re- 
coil (see Mic ic et all [2006. for more details). While the 
IMF of Pop III stars and BH growth are still a matter 
of debate, the goal in this paper is simply to examine 
the differences between 2LPT and ZA-seeded volumes. 
A different IMF or growth prescription will change ab- 
solute BH number densities and masses, and more (or 
less) aggressive BH growth will result in a difference in 
the ionized volumes. Hence, we chose the simplest Pop 
III IMF and BH growth prescription to contrast the two 

5 Note that our particle mass does not allow us to resolve all 
Pop Ill-hosting halos at the highest redshifts 



2LPT and the First Black Holes 



3 



Age of the Universe [Gyrs] 
0.9 0.5 0.3 0.2 

pc 1 1 



o 10 

Oh 



10 J r 



CM 

pq 



10' 



, 2 




I I I I I I I I I I I I I 




10 15 
redshift 



Fig. 1. — Top: The co- moving BH mass density as a function 
of redshift for 2LPT and ZA for the three boxes. The lines 
represent the mean of the 2LPT and ZA BH mass density, 
while the shaded region shows the maximum and minimum 
of the 2LPT and ZA quantities. Note the variation between 
the three boxes from cosmic variance. Bottom: the ratio of 
the 2LPT and ZA BH masses as a function of redshift. 2LPT 
BHs are more massive than the corresponding ZA BHs for 
the majority of cosmic time. 
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Fig. 2. — The probability density of the BH mass function at 
z — 6 computed by combining all three boxes for the 2LPT 
and ZA initialized volumes. The most massive black holes in 
the simulation volumes, log Mbh > 5.4, are consistently more 
likely in the 2LPT volumes by roughly 40%. 

initialization technique^. 

Once the BHs are seeded and accreting, we estimate 
the ionizing photon luminosity from the growing BH, as 
well as its effect on the neutral gas in the neighborhood 
through radiative transfer calculations. We assume that 
the BH is fed at the Eddington rate by major mergers, 
and this presupposes a rich supply of cold gas in these 
halos - an assumption that has not been observationally 
verified at such a high redshift. In fact, massive first 
stars can expel most of the gas in the host halo, starv- 

6 We assume that the radiative feedback from the accreting 
BH affects its growth by stopping the gas inflow after a Salpeter 
timescale 



ing the remnant BH dKitavama et~a l. 2004; Wh alen et al.l 
120041 : 1 Alvarez et a h 2009) until a merger with a gas-rich 
DMH can restock it. In particular, early simulations with 
~ 1OOM Pop III stars found that gas expulsion delayed 
accretion for 10 -50 Myr while the host DMH restocked 
its gas reservoir flJohnson et alJl2QQ7t IWise &; Abell2008l : 
I Wise et aT]|2012[ ). Recent simulations predict lower Pop 
III stellar masses, and their radiative feedback can still 
drive the majority of gas out of their host halos. How- 
ever, most of the gas infall at high redshift s and in rare 
peaks happen through cold streams - these cold streams 
present a small solid-angle to the stellar photon flux, 
which should make it more difficult to photo-evaporate. 
Thus, it may be plausible for Pop III DMHs, especially in 
more massive ones suppressed by LW radiation, to retain 
a reservoir of gas for future consumption. 

Our simulation is collisionless, so we calculate the im- 
pact of enhanced BH growth on the neutral gas numer- 
ically in post-processing. For an accreting BH, the tem- 
perature profile of the accretion disk is: 



T(r)=T (^) 3/4 (l-^ 



(6) 



/ 3GM BU M q 
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where, Mbh is BH mass, M gas is gas accretion rate, a 
is Stefan-Boltzmann's constant, r^ n is the innermost ac- 
cretion disk radius. The accretion disk luminosity at a 
given frequency, z/, can then be computed: 



2his 3 



exp 



hi 



kT(r) 



1 



rdr (7) 



where L v is the emitted luminosity per frequency band, 
c is the speed of light, k is Boltzmann's constant, and 
h is Planck's constant. The accretion efficiency, 77 con- 
verts the accreted mass (M acc ) to energy: M acc c 2 = 
l/2r g /r ini where r g = GMbu/c 2 . For a Schwarzschild 
BH r in is 



6r g with r] 



In practice, we used 



77 = 10% and adjusted r^ n accordingly. 

To quantify the effect of these BHs on reioniza- 
tion, we condu ct radiative transfe r calculations with 
Enzo+Moray ( Wise fc Abel 2011). After every snap- 
shot, we accumulate a fraction Q^/Qm of the mass into 
a fixed 512 3 grid, assuming all gas is hydrogen. At these 
scales and resolution, the gas should follow the DM. The 
density distribution is fixed during the radiative trans- 
fer calculation. We treat each accreting BH as a point 
source, with luminosities described by Equation (jfj) di- 
vided into bins of 13.6 - 40, 40 - 100, 100 - 10 3 , 10 3 - 10 4 
eV. Photon packages have the average photon energy in 
each bin, ~ 28, 70, 400, and 1200 eV, and they obey pe- 
riodic boundary c onditions. We use a nonequilibrium 
chemistry model (jAnninos et al.1 119971 ) with hydrogen 
only. We consider Compton cooling and free electron 
heating by the CMB, radiative losses from atomic cool- 
ing in the optically thin limit, and secondary ionizations 
from high-energy radiation, using the fitting formulae 
from lShull van Steenberd (|1985f ). 

3. RESULTS 
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Fig. 3. — Evolution of the mass- weighted ionization fraction 
of the ZA (solid) and 2LPT (dashed) simulations. The 2LPT- 
initialized simulation generates 20% more Pop Ill-hosting ha- 
los, which enhances BH merger rates and thus their total lu- 
minosity and contribution to reionization. 

By z = 6, we find that the 2LPT-initialized volume pro- 
duces 25more Pop III stars than the identical ZA volume. 
This is synonymous with saying that high mass DMHs 
collapsed earlier and more often. Indeed, at z = 24, 
there are 3 times the number of DMHs above 1O 7 M0, 
roughly a 5-cr peak, and by z = 20, the difference in the 
halo mass function above ~ 1O 7 M0 (~ 4 — a) is ~ 40%. 
Note that at z = 6, the difference in the number of newly 
formed Pop III star hosting halos has virtually vanished - 
this suggests that the simulation has undergone enough 
e- folding times to damp away transients. 

3.1. The Changes in the Seed Black Hole Population 

An overall 25% increase in the number of Pop Ill- 
host ing halos may not sound huge, but the consequences 
of this early and generous sowing of Pop III stars on 
proto-supermassive BH growth is more profound. When 
we grow seed BHs with the major-merger driven pre- 
scription outlined above, the minor difference in the halo 
mass function magnifies. Figure Q] shows the BH mass 
density, pbh, as a function of redshift. We calculate 
a factor of ~ 2 difference in pbh by z = 6. If we al- 
low a very vigorous, and perhaps unrealistic, BH growth 
scheme that fuels the BH continuously at the Eddington 
rate, then Pbh,2LPT would be more than 10 times larger 
than the ZA approximation. Future space-based gravita- 
tional wave observatories may detect these BH inspirals; 
the gravitational wave signal in this pre-reionization era 
could feature more numerous and louder sources than 
predicted, depending on the detector configuration. 

3.2. Implications for the Reionization Epoch 

Perhaps the most interesting consequence of the 2LPT 
DMH mass function is the effect of the more numerous 
and more massive BHs on reionization. Figure [2] shows 
the probability density of the BH mass function at z = 6 
over all 2LPT and ZA initialized volumes. We find that 
massive black holes (log Mbh ^ 5.4) are ~ 1.4 times 
more common in 2LPT volumes. 

Although our simulation is purely collisionless, we can 
estimate the contribution from the first BHs on reioniza- 




Fig. 4. — The response of neutral hydrogen in the Universe 
to the radiation produced by seed BHs in ZA (left) and 2LPT 
(right) initialized runs at three times in one realization. Both 
simulations used the same BH seeding and growth mecha- 
nism, and sampled the same initial phases. We used ray- 
tracing to propagate the BH radiation throughout the box, 
calculating the electron temperature of the gas and ionization 
state. 

tion. By z = 6, we find that 2LPT BHs emit 75% more 
ionizing photons per hydrogen atom than ZA BHs. This 
alone implies that the ionization history will be dramat- 
ically different. We next present the results from three- 
dimensional radiative transfer calculations on these six 
volumes. 

Figure [3] shows the evolution of the total BH luminos- 
ity and mass-weighted electron fraction in the ZA and 
2LPT simulations in one realization. At z > 10, the 
total luminosity fluctuates but is roughly equivalent. Af- 
terwards, the total luminosity in the 2LPT simulation 
is consistently ~ 3 times greater than ZA. This results 
in a larger ionized fraction at z < 9, and yields a to- 
tal ionized fraction of 0.16 and 0.25 for ZA and 2LPT at 
z = 6, respectively. However, this enhanced contribution 
to reionization only translates to a change in the Thom- 
son optical depth of Ar e = 1.7 x 10 -3 , if the universe is 
completely ionized at z < 6. In the other two realiza- 
tions, the ionized fraction increases by 0.12 and 0.09 at 
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z = 6. 

The density-weighted projections of electron fraction 
in Figure [4] depicts the spatial differences between ZA 
and 2LPT at three redshifts. The H II regions around 
the BHs are generally larger and more abundant in the 
2LPT case. Overall, the 2LPT ionized volume is ~ 2.5 
times larger, where we define ionized as x e > 0.5. Fur- 
thermore, this additional luminosity globally warms the 
IGM outside these H II to 17,000 K in underdense regions 
(S < 1), compared to 14,000 K in ZA. This IGM global 
warming by pre-reionization BHs may be important in 
regul ating the formation and growth of later, lower mass 
BHs (jTanaka et al.ll2012l ). 

4. DISCUSSION 

We found that pre-reionization simulations are partic- 
ularly sensitive to the phase-space initialization method, 
given that the astrophysically interesting DMHs- ones 
that host the first stars, BHs, and protogalaxies - are all 
expected to collapse at very high redshifts. This means 
that the transient errors in the initial positions and ve- 
locities have only a few e-folding times to damp away 
before we need to make reliable measurements of the 
collapsed DMHs. We found that the Zel'dovich approxi- 
mation underestimated the number and mass of the high 
mass DMHs during the dark ages, and this ~ 25% differ- 
ence in the extreme end of the halo mass function mag- 
nified, increasing the BH mass density ~ 1.4 times for 
log M B h > 5.4 at z = 6 (see Fig. [2j). This difference de- 
pends both on the fact that seed BHs form earlier, which 
allows them to grow for a longer time through mergers 
and accretion, and also that they are sown more often. 

One of the largest effects of this more accurate 2LPT 
technique is apparent in the reionization history of the 
early Universe. The volume of the Universe reionized at 
z = 6 more than doubles, only due to the initialization 
method. We caution that our main goal was to study 
differences in our two testbed simulations using simple 
prescriptions for BH feedback and mass accretion. We 
neglect cold gas inflows, mechanical feedback, and in- 
deed all gas and reionization physics is treated in post- 



processing. We are also neglecting a potentially critical 
effect of dark matter streaming motions on baryons at 
z = 1000, which will also delay galaxy formation com- 
pared to a standard approach (jTseliakhovich fc Hiratal 
l2Q10h . While the actual reionized volume of universe 
is not robust here, the critical point is that both ZA 
and Press- Schechter approaches underestimate the vol- 
ume merely because they both assume a linear evolution 
of DM over densities. 

We focused on Pop Ill-star hosting halos here, simply 
because the predicted DMH mass thresholds are easily 
defined. However, a ZA-initialized simulation will always 
underestimate the very high end of the DMH mass func- 
tion during the pre-reionization era, and this will have 
an effect on any prediction that relies on accurate high 
mass DMH number counts - an o bvious example is seed 
BH formation from direct collapse ( Bromm & Loeb 2003; 
lLodato fc Nataraiaril2^ iBegelman et alJl2QQd T 

Note, too, that our simulation is ~ 10 6 times too small 
to model billion solar mass BHs, whose number density 
is 1 per Gpc 3 ([Jiang et al.l l2009f ). Although WMAP-5 
can constrain the number of > 1O 13 M halos, our vol- 
ume is too small to sample these halos. We do pre- 
dict that these even rarer peaks will be more numer- 
ous and will be, on average, more massive in the pre- 
reionization epoch than predicted by N-body simulations 
initialized with ZA. In principle, observations will be able 
to constrain the true halo mass function in the dark 
ages. While the relatively small difference in r e we ob- 
served m ay not be disting uishable with PLANCK and 
Herschel (jZahn et al.l 120111 ). the g lobal IGM tempera - 
ture difference may be detectable (jTheuns et al.l [2QQ2h . 
and SKA could easily probe this d ifference in the 21cm 
power spectrum (Bar kana fc Loebll2005l ). 
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